library(tidyverse)
library(rio)
library(cregg)
library(patchwork)

# set replication folder as working directory
setwd("~replication")

load("data_genderedcost_long.rdata")

# only include completed answers - and answers given before deadline
# 2021-12-20 21:49:52 was the last response within the time frame
df_long <- df_long %>% 
  filter(SurveyStatus==2)

df_long <- df_long %>% 
  filter(SurveyEndTime<="2021-12-20 21:49:52")


####### FIGURE 4 ######
### MARGINAL MEANS FOR SUBSETS OF VICTIMS AND NON-VICTIMS

df_victims <- df_long %>% 
  dplyr::filter(victim=="Victim")

mm_victims <- cregg::mm(df_victims, # data
                        choice~position + remuneration + workload + work_environment, # formula
                        id = id)

mm_victims$victim <- "Victim" # add victim variable

df_nonvictims <- df_long %>% 
  dplyr::filter(victim=="Non-victim")

mm_nonvictims <- cregg::mm(df_nonvictims, # data
                           choice~position + remuneration + workload + work_environment, # formula
                           id = id)

mm_nonvictims$victim <- "Non-victim" # add victim variable

#combine estimates from the two subsets, add attribute variable for color and rearrange order of levels
mm_victim_subsets <- bind_rows(mm_victims, mm_nonvictims) %>% 
  mutate(feature = case_when(feature=="work_environment"~"Working Environment",
                             feature=="workload"~"Workload",
                             feature=="remuneration"~"Remuneration",
                             feature=="position"~"Position")) %>% 
  mutate(feature = factor(feature, levels = c("Position", "Remuneration", "Workload", "Working Environment")))

## LEFT HAND SIDE OF FIGURE 4
mm_victim_env_plot <- mm_victim_subsets %>% 
  filter(feature=="Working Environment") %>% 
  ggplot(data=., aes(y = level, x = estimate, color = victim, shape = victim)) +
  geom_point(position = position_dodge2(width = 0.5)) +
  geom_linerange(aes(xmin=estimate-std.error*1.39,
                     xmax=estimate+std.error*1.39),
                 position = position_dodge2(width = 0.5)) +
  xlab("Marginal mean") +
  ylab("") +
  scale_x_continuous(labels = seq(0.3,0.8,0.05), breaks = seq(0.3,0.8,0.05)) +
  geom_vline(xintercept = 0.5, linetype = "dashed") +
  theme_bw() +
  scale_color_grey("", start = 0.1, end = 0.65) +
  scale_shape_discrete("") +
  facet_wrap(~feature, scales = "free_y", ncol = 1) +
  theme(legend.position = "top", panel.background = element_rect(fill = "white"),
        strip.background = element_rect("white"),
        strip.text = element_text(hjust = 0, face = "bold"),
        panel.grid.major = element_blank(), panel.grid.minor = element_blank())



## estimate differences in MM for victims and nonvictims
mm_victim_differences <- cregg::mm_diffs(df_long, # data
                                         choice~position + remuneration + workload + work_environment, # formula
                                         ~victim)

mm_victim_differences <- mm_victim_differences %>% 
  mutate(feature = case_when(feature=="work_environment"~"Working Environment",
                             feature=="workload"~"Workload",
                             feature=="remuneration"~"Remuneration",
                             feature=="position"~"Position")) %>% 
  mutate(feature = factor(feature, levels = c("Position", "Remuneration", "Workload", "Working Environment")))


## RIGHT HAND SIDE IN FIGURE 4
mm_victim_env_dif_plot <- mm_victim_differences %>%
  filter(feature=="Working Environment") %>% 
  ggplot(data=., aes(y = level, x = estimate)) +
  geom_point(shape = 16, color = "black") +
  geom_linerange(aes(xmin=lower,
                     xmax=upper), color = "black") +
  xlab("Differences\n(victims - non-victims)") +
  ylab("") +
  scale_x_continuous(labels = seq(-0.2,0.2,0.05), breaks = seq(-0.2,0.2,0.05)) +
  geom_vline(xintercept = 0, linetype = "dashed") +
  theme_bw() +
  scale_color_grey("", start = 0.1, end = 0.65) +
  scale_shape_discrete("") +
  facet_wrap(~feature, scales = "free_y", ncol = 1) +
  coord_cartesian(xlim = c(-0.1,0.1)) +
  theme(legend.position = "top", panel.background = element_rect(fill = "white"),
        strip.background = element_rect("white"),
        strip.text = element_text(hjust = 0, face = "bold"),
        axis.text.y=element_blank(), axis.ticks.y=element_blank(),
        panel.grid.major = element_blank(), panel.grid.minor = element_blank())

## save in one plot
mm_victim_env_plot + mm_victim_env_dif_plot + plot_layout(widths = c(3, 1))

ggsave("figure4.pdf", height = 2.5, width = 10)
ggsave("figure4.png", height = 2.5, width = 10, dpi = 600)
